# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2020-2023, GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.

import os
import abc
import csv
import logging
import numpy
from shapely import geometry, wkt

from openquake.baselib.onnx import PicklableInferenceSession
from openquake.hazardlib import geo, imt, valid, InvalidFile
from openquake.sep.landslide.static_safety_factor import infinite_slope_fs
from openquake.sep.landslide.displacement import (
    critical_accel,
    jibson_2007_model_a,
    jibson_2007_model_b,
    cho_rathje_2022,
    fotopoulou_pitilakis_2015_model_a,
    fotopoulou_pitilakis_2015_model_b,
    fotopoulou_pitilakis_2015_model_c,
    fotopoulou_pitilakis_2015_model_d,
    saygili_rathje_2008,
    rathje_saygili_2009,
    jibson_etal_2000,
)
from openquake.sep.landslide.probability import(
    nowicki_jessee_2018,
    LANDCOVER_TABLE,
    LITHOLOGY_TABLE,
    allstadt_etal_2022_b,
    jibson_etal_2000_probability,
)
from openquake.sep.liquefaction.liquefaction import (
    hazus_liquefaction_probability,
    zhu_etal_2015_general,
    zhu_etal_2017_coastal,
    zhu_etal_2017_general,
    rashidian_baise_2020,
    allstadt_etal_2022,
    akhlagi_etal_2021_model_a,
    akhlagi_etal_2021_model_b,
    bozzoni_etal_2021_europe,
    todorovic_silva_2022_nonparametric_general,
    HAZUS_LIQUEFACTION_PGA_THRESHOLD_TABLE,
)
from openquake.sep.liquefaction.lateral_spreading import (
    hazus_lateral_spreading_displacement,
)
from openquake.sep.liquefaction.vertical_settlement import (
    hazus_vertical_settlement,
)
from os import path
import gzip


class SecondaryPeril(metaclass=abc.ABCMeta):
    """
    Abstract base class. Subclasses of SecondaryPeril have:

    1. a ``__init__`` method with global parameters coming from the job.ini
    2. a ``prepare(sitecol)`` method modifying on the site collection, called
    in the ``pre_execute`` phase, i.e. before running the calculation
    3. a ``compute(mag, imt, gmf, sites)`` method called during the calculation
    of the GMFs; gmf is an array of length N1 and sites is a (filtered)
    site collection of length N1 (normally N1 < N, the total number of sites)
    4. an ``outputs`` attribute which is a list of column names which will be
    added to the gmf_data array generated by the ground motion calculator

    The ``compute`` method will return a tuple with ``O`` arrays where ``O``
    is the number of outputs.
    """
    peril = ''  # overriden in subclasses
    outputs = []

    @classmethod
    def __init_subclass__(cls):
        # make sure the name of the outputs are valid IMTs
        for out in cls.outputs:
            imt.from_string(out)

    @classmethod
    def instantiate(cls, secondary_perils, sec_peril_params, oq):
        """
        :param secondary_perils: a list of class names
        :param sec_peril_params: a list of dicts with instantiation params
        :param oq: OqParam object to be attached to the instances
        :returns: a list of instances of the secondary peril classes
        """
        if not sec_peril_params:
            sec_peril_params = [{}] * len(secondary_perils)
        instances = []
        for clsname, params in zip(secondary_perils, sec_peril_params):
            obj = globals()[clsname](**params)
            obj.oq = oq
            instances.append(obj)
        return instances

    @abc.abstractmethod
    def prepare(self, sites):
        """Add attributes to sites"""

    @abc.abstractmethod
    def compute(self, mag, imt_gmf, sites):
        """
        :param mag: magnitude
        :param imt_gmf: a list of pairs (imt, gmf)
        :param sites: a filtered site collection
        """

    def __repr__(self):
        return "<%s>" % self.__class__.__name__


class HazusLiquefaction(SecondaryPeril):
    peril = 'liquefaction'
    outputs = ["LiqProb"]

    def __init__(self, map_proportion_flag=True):
        self.map_proportion_flag = map_proportion_flag

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                out.append(
                    hazus_liquefaction_probability(
                        pga=gmf,
                        mag=mag,
                        liq_susc_cat=sites.liq_susc_cat,
                        groundwater_depth=sites.gwd,
                        do_map_proportion_correction=self.map_proportion_flag))
        return out


class HazusDeformation(SecondaryPeril):
    """
    Computes PGDMax or PGDGeomMean from PGA
    """
    peril = 'liquefaction'
    outputs = ["PGDMax"]

    def __init__(
        self,
        return_unit="m",
        deformation_component="PGDMax",
        pga_threshold_table=HAZUS_LIQUEFACTION_PGA_THRESHOLD_TABLE,
    ):
        self.return_unit = return_unit
        self.deformation_component = getattr(imt, deformation_component)
        self.outputs = [deformation_component]

        if pga_threshold_table != HAZUS_LIQUEFACTION_PGA_THRESHOLD_TABLE:
            pga_threshold_table = {
                k.encode("utf-8"): v for k, v in pga_threshold_table.items()
            }
        self.pga_threshold_table = pga_threshold_table

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                ls = hazus_lateral_spreading_displacement(
                    mag=mag,
                    pga=gmf,
                    liq_susc_cat=sites.liq_susc_cat,
                    pga_threshold_table=self.pga_threshold_table,
                    return_unit=self.return_unit,
                )
                vs = hazus_vertical_settlement(
                    sites.liq_susc_cat, return_unit=self.return_unit
                )
                out.append(self.deformation_component(ls, vs))
        return out


class ZhuEtAl2015LiquefactionGeneral(SecondaryPeril):
    """
    Computes the liquefaction probability from PGA and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur"]

    def __init__(
        self,
        intercept=24.1,
        pgam_coeff=2.067,
        cti_coeff=0.355,
        vs30_coeff=-4.784,
    ):
        self.intercept = intercept
        self.pgam_coeff = pgam_coeff
        self.cti_coeff = cti_coeff
        self.vs30_coeff = vs30_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                prob_liq, out_class = zhu_etal_2015_general(
                    pga=gmf, mag=mag, cti=sites.cti, vs30=sites.vs30)
                out.append(prob_liq)
                out.append(out_class)
        return out


class ZhuEtAl2017LiquefactionCoastal(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur", "LSE"]

    def __init__(
        self,
        intercept=12.435,
        pgv_coeff=0.301,
        vs30_coeff=-2.615,
        dr_coeff=0.0666,
        dc_coeff=-0.0287,
        dcdr_coeff=-0.0369,
        precip_coeff=0.0005556,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.vs30_coeff = vs30_coeff
        self.dr_coeff = dr_coeff
        self.dc_coeff = dc_coeff
        self.dcdr_coeff = dcdr_coeff
        self.precip_coeff = precip_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                prob_liq, out_class, lse = zhu_etal_2017_coastal(
                    pgv=gmf,
                    vs30=sites.vs30,
                    dr=sites.dr,
                    dc=sites.dc,
                    precip=sites.precip)
                out.append(prob_liq)
                out.append(out_class)
                out.append(lse)
        return out


class ZhuEtAl2017LiquefactionGeneral(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur", "LSE"]

    def __init__(
        self,
        intercept=8.801,
        pgv_scaling_factor=1.0,
        pgv_coeff=0.334,
        vs30_coeff=-1.918,
        dw_coeff=-0.2054,
        wtd_coeff=-0.0333,
        precip_coeff=0.0005408,
    ):
        self.intercept = intercept
        self.pgv_scaling_factor = pgv_scaling_factor
        self.pgv_coeff = pgv_coeff
        self.vs30_coeff = vs30_coeff
        self.dw_coeff = dw_coeff
        self.wtd_coeff = wtd_coeff
        self.precip_coeff = precip_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                prob_liq, out_class, lse = zhu_etal_2017_general(
                    pgv=gmf,
                    vs30=sites.vs30,
                    dw=sites.dw,
                    wtd=sites.gwd,
                    precip=sites.precip)
                out.append(prob_liq)
                out.append(out_class)
                out.append(lse)
        return out


class RashidianBaise2020Liquefaction(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and PGA and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur", "LSE"]

    def __init__(
        self,
        intercept=8.801,
        pgv_scaling_factor=1.0,
        pgv_coeff=0.334,
        vs30_coeff=-1.918,
        dw_coeff=-0.2054,
        wtd_coeff=-0.0333,
        precip_coeff=0.0005408,
    ):
        self.intercept = intercept
        self.pgv_scaling_factor = pgv_scaling_factor
        self.pgv_coeff = pgv_coeff
        self.vs30_coeff = vs30_coeff
        self.dw_coeff = dw_coeff
        self.wtd_coeff = wtd_coeff
        self.precip_coeff = precip_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute liquefaction "
                "probability using the RashidianBaise2020Liquefaction model")
        prob_liq, out_class, lse = rashidian_baise_2020(
            pga=pga,
            pgv=pgv,
            vs30=sites.vs30,
            dw=sites.dw,
            wtd=sites.gwd,
            precip=sites.precip,
        )
        out.append(prob_liq)
        out.append(out_class)
        out.append(lse)
        return out


class AllstadtEtAl2022Liquefaction(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and PGA and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur", "LSE"]

    def __init__(
        self,
        intercept=8.801,
        pgv_coeff=0.334,
        vs30_coeff=-1.918,
        dw_coeff=-0.2054,
        wtd_coeff=-0.0333,
        precip_coeff=0.0005408,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.vs30_coeff = vs30_coeff
        self.dw_coeff = dw_coeff
        self.wtd_coeff = wtd_coeff
        self.precip_coeff = precip_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute liquefaction "
                "probability using the AllstadtEtAl2022Liquefaction model"
            )

        prob_liq, out_class, lse = allstadt_etal_2022(
            pga=pga,
            pgv=pgv,
            mag=mag,
            vs30=sites.vs30,
            dw=sites.dw,
            wtd=sites.gwd,
            precip=sites.precip,
        )
        out.append(prob_liq)
        out.append(out_class)
        out.append(lse)
        return out


class AkhlagiEtAl2021LiquefactionA(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and transforms it
    to binary output via the predefined probability threshold.
    """
    experimental = True
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur"]

    def __init__(
        self,
        intercept=4.925,
        pgv_coeff=0.694,
        tri_coeff=-0.459,
        dc_coeff=-0.403,
        dr_coeff=-0.309,
        zwb_coeff=-0.164,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.tri_coeff = tri_coeff
        self.dc_coeff = dc_coeff
        self.dr_coeff = dr_coeff
        self.zwb_coeff = zwb_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                prob_liq, out_class = akhlagi_etal_2021_model_a(
                    pgv=gmf,
                    tri=sites.tri,
                    dc=sites.dc,
                    dr=sites.dr,
                    zwb=sites.zwb)
                out.append(prob_liq)
                out.append(out_class)
        return out


class AkhlagiEtAl2021LiquefactionB(SecondaryPeril):
    """
    Computes the liquefaction probability from PGV and transforms it
    to binary output via the predefined probability threshold.
    """
    experimental = True
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur"]

    def __init__(
        self,
        intercept=9.504,
        pgv_coeff=0.706,
        vs30_coeff=-0.994,
        dc_coeff=-0.389,
        dr_coeff=-0.291,
        zwb_coeff=-0.205,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.vs30_coeff = vs30_coeff
        self.dc_coeff = dc_coeff
        self.dr_coeff = dr_coeff
        self.zwb_coeff = zwb_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                prob_liq, out_class = akhlagi_etal_2021_model_b(
                    pgv=gmf,
                    vs30=sites.vs30,
                    dc=sites.dc,
                    dr=sites.dr,
                    zwb=sites.zwb)
                out.append(prob_liq)
                out.append(out_class)
        return out


class Bozzoni2021LiquefactionEurope(SecondaryPeril):
    """
    Computes the liquefaction probability from PGA and transforms it
    to binary output via the predefined probability threshold.
    """
    peril = 'liquefaction'
    outputs = ["LiqProb", "LiqOccur"]

    def __init__(
        self,
        intercept=-11.489,
        pgam_coeff=3.864,
        cti_coeff=2.328,
        vs30_coeff=-0.091,
    ):
        self.intercept = intercept
        self.pgam_coeff = pgam_coeff
        self.cti_coeff = cti_coeff
        self.vs30_coeff = vs30_coeff

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                prob_liq, out_class = bozzoni_etal_2021_europe(
                    pga=gmf, mag=mag, cti=sites.cti, vs30=sites.vs30)
                out.append(prob_liq)
                out.append(out_class)
        return out


class TodorovicSilva2022NonParametric(SecondaryPeril):
    """
    Computes the liquefaction occurrence from PGV and transforms it
    to probability of belonging to positive class, i.e., liquefaction
    occurrence.
    """
    peril = 'liquefaction'
    outputs = ["LiqOccur", "LiqProb"]

    def prepare(self, sites):
        model_file = "liquefaction/data/todorovic_silva_2022/" + \
            "todorovic_silva_2022.onnx.gz"
        model_path = path.join(path.dirname(__file__), model_file)
        with gzip.open(model_path, "rb") as f:
            self.model = f.read()
        self.inference_session = PicklableInferenceSession(self.model)

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                out_class, out_prob = \
                    todorovic_silva_2022_nonparametric_general(
                        pgv=gmf,
                        vs30=sites.vs30,
                        dw=sites.dw,
                        wtd=sites.gwd,
                        precip=sites.precip,
                        session=self.inference_session)
                out.append(out_class)
                out.append(out_prob)
        return out


supported = [cls.__name__ for cls in SecondaryPeril.__subclasses__()]


class Jibson2007ALandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides
    as function of pga and critical acceleration.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def __init__(
            self, c1=0.215, c2=2.341, c3=-1.438, crit_accel_threshold=0.05):
        self.c1 = c1
        self.c2 = c2
        self.c3 = c3
        self.crit_accel_threshold = crit_accel_threshold

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                Disp = jibson_2007_model_a(
                    gmf, sites.crit_accel,
                    self.c1, self.c2, self.c3,
                    self.crit_accel_threshold)
                out.append(Disp)
        return out


class Jibson2007BLandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides
    as function of magnitude, pga and critical acceleration.
    It is recommended when magnitude is betweeen 
    5.3 and 7.6.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def __init__(
        self,
        c1=-2.71,
        c2=2.335,
        c3=-1.478,
        c4=0.424,
        crit_accel_threshold=0.05,
    ):
        self.c1 = c1
        self.c2 = c2
        self.c3 = c3
        self.c4 = c4
        self.crit_accel_threshold = crit_accel_threshold

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                Disp = jibson_2007_model_b(
                    gmf, sites.crit_accel, mag,
                    self.c1, self.c2, self.c3, self.c4,
                    self.crit_accel_threshold)
                out.append(Disp)
        return out


class ChoRathje2022Landslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides from
    pgv and considering the slope fundamental period.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                Disp = cho_rathje_2022(
                        gmf,
                        sites.tslope,
                        sites.crit_accel,
                        sites.hratio)
                out.append(Disp)
        return out


class FotopoulouPitilakis2015ALandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides from
    pgv and moment magnitude.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                Disp = fotopoulou_pitilakis_2015_model_a(
                    gmf, mag, sites.crit_accel)
                out.append(Disp)
        return out      
        
        
class FotopoulouPitilakis2015BLandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides from
    pga and moment magnitude.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                Disp = fotopoulou_pitilakis_2015_model_b(
                          gmf, mag, sites.crit_accel,)
                out.append(Disp)
        return out      
        
class FotopoulouPitilakis2015CLandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements from pga (in terms of ratio with
    the landslide critical acceleration) and moment magnitude.
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                Disp = fotopoulou_pitilakis_2015_model_c(
                    gmf, mag, sites.crit_accel)
                out.append(Disp)
        return out


class FotopoulouPitilakis2015DLandslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements of landslides from pgv and 
    pga
    '''    
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute landslide disp "
                "according to Fotopoulou_Pitilakis_2015_PGV_PGA")

        Disp = fotopoulou_pitilakis_2015_model_d(
            pgv,
            pga,
            sites.crit_accel,
            )
        out.append(Disp)
        return out 


class SaygiliRathje2008Landslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements from pga and pgv
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope))
        
    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute landslide "
                "disp according to Saygili_Rathje_2008")

        Disp = saygili_rathje_2008(pga, pgv, sites.crit_accel)
        out.append(Disp)
        return out


class RathjeSaygili2009Landslides(SecondaryPeril):   
    '''
    Computes earthquake-induced displacements from pga and moment 
    magnitude
    '''
    peril = 'landslide'
    outputs = ["Disp"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope))
        
        print(sites)

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "PGA":
                Disp = rathje_saygili_2009(gmf, mag, sites.crit_accel)
                out.append(Disp)
        return out  


class JibsonEtAl2000Landslides(SecondaryPeril):
    '''
    Computes earthquake-induced displacements and related probability according
    to Jibson et al. (2000) as function of arias intensity.
    '''
    peril = 'landslide'
    outputs = ["Disp", "DispProb"]

    def prepare(self, sites):
        sites.add_col(
            "Fs",
            float,
            infinite_slope_fs(
                slope=sites.slope,
                cohesion=sites.cohesion_mid,
                friction_angle=sites.friction_mid,
                saturation_coeff=sites.saturation,
                soil_dry_density=sites.dry_density,
                slab_thickness=sites.slab_thickness,
            ),
        )
        sites.add_col(
            "crit_accel", float, critical_accel(sites.Fs, sites.slope)
        )

    def compute(self, mag, imt_gmf, sites):
        out = []
        for im, gmf in imt_gmf:
            if im.string == "IA":
                Disp = jibson_etal_2000(
                    gmf, sites.crit_accel)
                out.append(Disp)
                out.append(jibson_etal_2000_probability(Disp))
        return out        


class NowickiJessee2018Landslides(SecondaryPeril):
    """
    Computes the landslide probability from PGV and areal coverage.
    """
    peril = 'landslide'
    outputs = ["LsProb", "LSE"]

    def __init__(
        self,
        intercept: float = -6.30,
        pgv_coeff: float = 1.65,
        slope_coeff: float = 0.06,
        coeff_table_lith = LITHOLOGY_TABLE,
        coeff_table_cov = LANDCOVER_TABLE,
        cti_coeff: float = 0.03,
        interaction_term: float = 0.01,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.slope_coeff = slope_coeff
        self.coeff_table_lith = coeff_table_lith
        self.coeff_table_cov = coeff_table_cov
        self.cti_coeff = cti_coeff
        self.interaction_term = interaction_term

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute landslide "
                "probability using the NowickiJessee2018Landslides model"
            )
        
        prob_ls, lse = nowicki_jessee_2018(
            pgv = pgv,
            slope=sites.slope,
            lithology=sites.lithology,
            landcover=sites.landcover,
            cti=sites.cti,
        )
        out.append(prob_ls)
        out.append(lse)
            
        return out


class AllstadtEtAl2022Landslides(SecondaryPeril):
    """
    Corrects LSE according to Allstadt et al. (2022).
    """
    peril = 'landslide'
    outputs = ["LsProb", "LSE"]

    def __init__(
        self,
        intercept: float = -6.30,
        pgv_coeff: float = 1.65,
        slope_coeff: float = 0.06,
        coeff_table_lith = LITHOLOGY_TABLE,
        coeff_table_cov = LANDCOVER_TABLE,
        cti_coeff: float = 0.03,
        interaction_term: float = 0.01,
    ):
        self.intercept = intercept
        self.pgv_coeff = pgv_coeff
        self.slope_coeff = slope_coeff
        self.coeff_table_lith = coeff_table_lith.copy()
        self.coeff_table_lith["su"] = -1.36
        self.coeff_table_lith[b"su"] = -1.36 
        self.coeff_table_cov = coeff_table_cov
        self.cti_coeff = cti_coeff
        self.interaction_term = interaction_term

    def prepare(self, sites):
        pass

    def compute(self, mag, imt_gmf, sites):
        out = []
        pga = None
        pgv = None
        for im, gmf in imt_gmf:
            if im.string == "PGV":
                pgv = gmf
            elif im.string == "PGA":
                pga = gmf
            else:
                continue
        # Raise error if either PGA or PGV is missing
        if pga is None or pgv is None:
            raise ValueError(
                "Both PGA and PGV are required to compute landslide "
                "probability using the AllstadtEtAl2022Landslides model"
            )
        
        prob_ls, lse = allstadt_etal_2022_b(
            pga = pga,
            pgv = pgv,
            slope=sites.slope,
            lithology=sites.lithology,
            landcover=sites.landcover,
            cti=sites.cti,
        )
        
        
        out.append(prob_ls)
        out.append(lse)
            
        return out


def csv2peril(fname, name, sitecol, tofloat, asset_hazard_distance):
    """
    Converts a CSV file into a peril array of length N
    """
    data = []
    with open(fname) as f:
        for row in csv.DictReader(f):
            intensity = tofloat(row['intensity'])
            if intensity > 0:
                data.append((valid.longitude(row['lon']),
                             valid.latitude(row['lat']),
                             intensity))
    data = numpy.array(data, [('lon', float), ('lat', float),
                              ('number', float)])
    logging.info('Read %s with %d rows' % (fname, len(data)))
    if len(data) != len(numpy.unique(data[['lon', 'lat']])):
        raise InvalidFile('There are duplicated points in %s' % fname)
    try:
        distance = asset_hazard_distance[name]
    except KeyError:
        distance = asset_hazard_distance['default']
    sites, filtdata, _discarded = geo.utils.assoc(
        data, sitecol, distance, 'filter')
    peril = numpy.zeros(len(sitecol), float)
    peril[sites.sids] = filtdata['number']
    return peril


def wkt2peril(fname, name, sitecol):
    """
    Converts a WKT file into a peril array of length N
    """
    with open(fname) as f:
        header = next(f)  # skip header
        if header != 'geom\n':
            raise ValueError('%s has header %r, should be geom instead' %
                             (fname, header))
        text = f.read()
        if not text.startswith('"'):
            raise ValueError('The geometry must be quoted in %s : "%s..."' %
                             (fname, text.split('(')[0]))
        geom = wkt.loads(text.strip('"\n'))  # strip quotes and newlines
    peril = numpy.zeros(len(sitecol), float)
    for sid, lon, lat in sitecol.complete.array[['sids', 'lon', 'lat']]:
        peril[sid] = geometry.Point(lon, lat).within(geom)
    return peril



class Volcanic(SecondaryPeril):
    """
    Import ASH, LAVA, LAHAR, PYRO from CSV files
    """
    peril = 'volcanic'
    outputs = ["ASH", "LAVA", "LAHAR", "PYRO"]

    def prepare(self, sites):
        """
        Import the CSV files for the volcanic subperils
        """
        for peril in self.oq.inputs['multi_peril']:
            assert peril in self.outputs, peril
        self.fname_by_peril = self.oq.inputs['multi_peril']
        N = len(sites)
        self.data = {'sid': sites.sids, 'eid': numpy.zeros(N, numpy.uint32)}
        names = []
        for name, fname in self.fname_by_peril.items():
            fname = os.path.join(self.oq.base_path, fname)
            tofloat = (valid.positivefloat if name == 'ASH'
                       else valid.probability)
            with open(fname) as f:
                header = next(f)
            if 'geom' in header:
                peril = wkt2peril(fname, name, sites)
            else:
                peril = csv2peril(fname, name, sites, tofloat,
                                  self.oq.asset_hazard_distance)
            if peril.sum() == 0:
                logging.warning('No sites were affected by %s' % name)
            self.data[f'{self.__class__.__name__}_{name}'] = peril
            names.append(name)

    def compute(self, mag, imt_gmf, sites):
        # doing nothing, since all the work is in the `prepare` method
        return []


LIQUEFACTION_MODELS = {cls.__name__ for cls in SecondaryPeril.__subclasses__()
                       if cls.peril == 'liquefaction'}
LANDSLIDE_MODELS = {cls.__name__ for cls in SecondaryPeril.__subclasses__()
                    if cls.peril == 'landslide'}


def corresponds(col, peril, imt):
    """
    Associate the given IMT to a column in gmf_data
    """
    if not col.endswith(imt):
        return False
    if peril == 'groundshaking':
        return True
    name, _ = col.split('_')
    if peril == 'liquefaction':
        return name in LIQUEFACTION_MODELS
    elif peril == 'landslide':
        return name in LANDSLIDE_MODELS
    else:
        raise NameError(f'Unknown peril {peril}')
